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Abstract. A sharp stability analysis of atomistic-to-continuum coupling methods 
is essential for evaluating their capabilities for predicting the formation and motion 
of lattice defects. We formulate a simple one-dimensional model problem and give a 
detailed analysis of the stability of the force-based quasicontinuum (QCF) method. 
The focus of the analysis is the question whether the QCF method is able to predict 
a critical load at which fracture occurs. 

Numerical experiments show that the spectrum of a linearized QCF operator is 
identical to the spectrum of a linearized energy-based quasi-nonlocal quasicontinuum 
operator (QNL), which we know from our previous analyses to be positive below the 
critical load. However, the QCF operator is non-normal and it turns out that it is 
not generally positive definite, even when all of its eigenvalues are positive. Using a 
combination of rigorous analysis and numerical experiments, we investigate in detail 
for which choices of "function spaces" the QCF operator is stable, uniformly in the 
size of the atomistic system. 

Force-based multi-physics coupling methods are popular techniques to circum- 
vent the difficulties faced in formulating consistent energy-based coupling approaches. 
Even though the QCF method is possibly the simplest coupling method of this kind, 
we anticipate that many of our observations apply more generally. 



1. Introduction 

Low energy equilibria for crystalline materials are typically characterized by local- 
ized defects that interact with their environment through long-range elastic fields. 
Atomistic-to-continuum coupling methods seek to make the accurate computation of 
such problems possible by using the accuracy of atomistic modeling only in the neigh- 
borhood of defects where the deformation is highly non-uniform. At some distance 
from the defects, sufficient accuracy can be obtained by the use of continuum models, 
which facilitate the reduction of degrees of freedom. The accuracy of the atomistic 
model at the defect combined with the efficiency of a continuum model for the far field 
enables, in principle, the reliable simulation of systems that are inaccessible to pure 
atomistic or pure continuum models. 
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Typical test problems for atomistic-to-continuum coupling methods have been dis- 
location formation under an indenter, crack tip deformation, and deformation and 
fracture of grain boundaries [18]. In each of these problems, the crystal deforms qua- 
sistatically until the equilibrium equations become singular, for example, when a dis- 
location is formed or moves or when a crack tip advances. Depending of the nature of 
the singularity, the crystal will then typically undergo a dynamic process when further 
loaded. 

The quasicontinuum (QC) approximation models the continuum region by using an 
energy density that exactly reproduces the lattice-based energy density at uniform 
strain (the Cauchy-Born rule). Several variants of the QC approximation have been 
proposed that differ in how the atomistic and continuum regions are coupled [3,5,11,13, 
14,18,20,24,26]. Analyses of QC approximation have been given in [7,12,16,17,19,21]. 
We refer to [8] for a detailed review of the formulation and analysis, relevant to the 
present work, of different QC methods. Other coupling models are analyzed in [1,22,23]. 

In [8], we have begun to investigate whether the QC method can reliably predict 
the formation of defects. The main ingredient to establish whether or not this is the 
case is a sharp analysis to predict under which conditions the QC method is "stable." 
More precisely, we ask whether there exist "stable" solutions of the QC method up 
to a critical load. We have begun to investigate this question in some depth for the 
most common energy-based QC formulations in [8]. In the present paper, we present 
a corresponding sharp stability analysis for the force-based quasicontinuum (QCF) 
method [4,5,24]. 

We focus on a one dimensional periodic chain with next-nearest neighbour pair inter- 
actions, which is introduced in Section 12.11 For this model, the uniform configuration 
ceases to be stable when the applied tensile strain reaches a critical value (fracture). 

For the atomistic model and for energy-based QC formulations, coercivity (positiv- 
ity) of the second variation evaluated at the equilibrium solution provides the natural 
notion of stability. However, the QCF method, which we describe in Sections 12.31 and 
12.51 leads to non-conservative equilibrium equations, and therefore, positivity of the 
linearized QCF operator may be an inappropriate notion of stability. Indeed, we prove 
in Section I4TT1 that, generically, the linearized QCF operator is indefinite. 

As a consequence, we consider two further notions of stability. First, we investigate 
for which choices of discrete function spaces (that is, for which choices of topologies) 
does the linearized QCF operator have an inverse that is bounded uniformly in the size 
of the atomistic system. In Section 14.21 we present several sharp stability results as 
well as interesting counterexamples. However, these operator stability results do not 
necessarily correspond to any physical notion of stability. Hence, in Section 14.41 we 
propose the notion of dynamical stability, which can be reduced to certain properties of 
the eigenvalues. A careful numerical study suggests that the spectrum of the linearized 
QCF operator and that of the linearized quasi-nonlocal QC operator (QNL) (see [26] 
and Section 14731 are identical. Combined with our previous results [8], this indicates 
that the QCF method is dynamically stable up to the critical load for fracture. 
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2. The force-based quasicontinuum method 

2.1. The atomistic model problem. We consider deformations from the reference 
lattice eZ, where e > is a scaling that we will fix below. For the sake of simplicity, 
we admit only deformations which are periodic displacements from the uniform state 
Up = FeL = (Fs£)e e z, that is, we admit deformations from the space 

yF = Uf where 

U = {u G K z : u e+2 N = u e for £ G Z, and J2f=-N+i u t = °}- 

We call F the macroscopic deformation gradient, and we set e = 1/N throughout. 
Although the energies and forces are defined for general 2iV-periodic displacements, 
we only admit those with zero mean, as is common for continuum problems with 
periodic boundary conditions, in order to obtain unique solutions to the equilibrium 
equations. 

We consider only nearest-neighbor and next-nearest neighbor pair interactions so 
that the potential energy per period of a deformation y G is given by 

TV 

£M = £ W$ + 0(^ + ^+i))> 

where 

Ve = ^(ye-Ve-i), 
and where is a Lennard- Jones type interaction potential: 

(i) 0GC 3 ((O,+oo);M), 

(ii) there exists r* > such that is convex in (0, r*) and concave in (r*, +oo). 

(iii) 4>( k \r) — > rapidly as r / oo, for k — 0, . . . , 3. 

Assumption (iii) is not strictly necessary for our analysis but serves to motivate that 
next-nearest neighbour interactions are typically dominated by nearest-neighbour terms. 

We assume that the atomistic system is subject to 2iV-periodic external forces (fejeez 
with zero mean, i.e., / G U, so that the total energy per period takes the form 

TV 

£ t a ° t {y)=£ a {y)-E £ f m . 

e=-N+i 

Equilibria y G of the atomistic total energy are solutions to the equilibrium 
equations 

F*s(y) + ft = 0, -oo < £ < oo, (1) 
where the (scaled) atomistic forces jF a : y F U* are defined by 

1 dS a (y) 
£ dye 

and where U* is the space of linear functionals on U. We remark that the translational 
invariance of the atomistic energy implies that J-^t{v) nas zero mean > 

N , 

e=-N+i 



4 



M. DOBSON, M. LUSKIN, AND C. ORTNER 



(3) 



(4) 



where e = {l)eez is the unit translation vector. Thus, we see that, at least heuristically, 
the external force vector / lies indeed in the range of the atomistic force operator. 
We note, moreover, that yp is an equilibrium of the atomistic energy, that is 

FhAvf) = - oo < £ < oo, for all F > 0. 

The question which we will investigate in this paper, beginning in Section [31 is for 
which F it is a stable equilibrium, and whether the force-based QC method is able to 
predict the stability of yp. 

2.2. The local quasicontinuum approximation. We begin by observing that the 
atomistic energy can be rewritten as a sum over the contributions from each atom, 

N 

S a (y) = e E i(v) where 

e=-N+i 

E e(y) = \ Wi) + <f>(y'i + i) + M-i + vi) + </>(y' t+ i + <4 +2 )] • 

If y is "smooth" , that is, if y[ varies slowly, then the atomistic energy can be accurately 
approximated by the Cauchy-Born or local quasicontinuum energy 

N 

S qc i(y) = £ E t(y)> where 

e=-N+i 

E c e (y) = § + <t>{y' i+x ) + <t>(2y' e ) + <f>{2j/ t+1 )] = \ [<j> ch {y'e) + <!>M+i)\ , 

where c b( r ) = 4'{ r ) + 0(2r) is the Cauchy-Born stored energy density. 

In this approximation we have replaced the next-nearest neighbor interactions by 
nearest neighbor interactions to obtain a model with stronger locality. This makes it 
possible to coarsen the model (to remove degrees of freedom), which eventually leads to 
significant gains in efficiency [5,18]. However, in the present work we will not consider 
this additional step. 

An equilibrium y e yp of the local QC energy is a solution to the equilibrium 
equations 

FcAv) + h = °> -oo < £ < oo, (5) 
where the (scaled) local QC forces JF C : y F — > U* are defined by 

1 d£ & {y) 

s C; t{y) := E , -oo < £ < oo. 

£ dye 

As in ([2]) it follows that the vector J- C {y) has zero mean. 

2.3. The force-based quasicontinuum approximation. If a deformation y is "smooth" 
except in a small region of the domain, then it is desirable to couple the accurate 
atomistic description with the efficient continuum description. The force-based qua- 
sicontinuum (QCF) approximation achieves this by mixing the equilibrium equations 

of the atomistic model with those of the continuum model without any interface or 
transition region. 

Suppose that y is "smooth" except in a region A := {—K, . . . ,K}, where K > 1. 
We call A the atomistic region and C = {—N + 1, . . . , N} \ A the continuum region. 
The force-based QC approximation is obtained by evaluating the forces in the atomistic 
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region by the full atomistic model ([T]) and the forces in the continuum region by the local 
QC model (jSJ). This yields the QCF operator for the (scaled) forces .F qc f : — * U* ■, 
defined by 

F^AV) - { jr cAyl if£GC . (6) 

Force-based coupling methods such as ([6]) are trivially consistent (provided the con- 
tinuum model is consistent with the atomistic model) and are therefore a natural rem- 
edy for the inconsistencies one observes when formulating simple energy-based coupling 
methods such as the original QC method [20]. Similar constructions have appeared in 
the literature under several different names and for various applications (e.g., FeAt [15], 
CADD [25], or brutal force mixing [2]). In the context of the QC method this method 
was first described in [5], where it was shown that the force-based QC method is the 
limit of the so-called ghost-force correction iteration [24] . A basin of attraction and rate 
for the convergence the ghost-force correction iteration to the force-based QC method 
was given in [5]. Sharp stability estimates for the ghost-force correction iteration are 
given in [10]. 

Unfortunately, the forces generated by the QCF method are non-conservative, and 
hence cannot be associated with an energy. Moreover, even though both the atomistic 
forces J- Ay) an d the local QC forces J- C {y) have zero mean, it turns out that this is 
false for the mixed forces ^ r qc f(?/). A straightforward computation shows that 

N 

£ r**M = e ~W(2yi.*) - W-s + y'-K-i) - W-K+i + y'-K)] 

e=-N+i 

- e- 1 [2(p\2y' K+1 ) - (f)'(y' K+2 + y' K+1 ) - <j)'{y' K+l + y' K )}, 

which is in general non-zero. After introducing the necessary notation, we will overcome 
this difficulty by defining a variational form of the QCF method, which effectively 
projects the QCF forces onto the correct range. 

2.4. Norms and variational notation. For future reference, we recall the backward 
first difference v' e = e~ 1 (f^ — i^-i) and also define the centered second difference v" t = 
£~ 2 {vt+i -2v e + v e _i). 

For displacements v G U and 1 < p < oo, we define the i p e norms, 




1 < p < oo, 



^ max^jv+i,...,^ 1^1, p = oo, 

and let U°' p denote the space U equipped with the if. norm. We further define the W 1,p 
norm 

Hli/ 1 * : = IKIk> 

and let W 1,p denote the space U equipped with the U x ' p norm. Similarly, we define the 
space U 2 ' p and its associated U 2,p norm. 

The inner product associated with the l 2 e norm is 

N 

(v,w):=e ViW£ for v,w eU. 

e=-N+i 
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We have defined the norms || • \\ £ p and the inner product (•, •) on U, though we will also 
apply them for arbitrary vectors from M. 2N . 

The external force / = {fe)e^z is a 2iV-periodic mean zero vector, and we have seen 
that the atomistic forces and the forces in the QCL method are also 2iV-periodic mean 
zero vectors. Using the inner product, we can view fg as a linear functional on U. We 
recall that the space of linear functionals on 14 is denoted by U*, and we note that each 
such TeW* has a unique representation as a zero mean 2iV-periodic vector G U, 

T[v] = (gr,v) Vt> G IA. (7) 

We will normally not make a distinction between these representations. For example, 
an external force vector / may be equally interpreted as a linear functional (i.e., / G 
U*), or identified with its Riesz representation (i.e., / G U). 

For g G U*, s = 0, 1, and 1 < p < oo, we define the negative norms ||p||w- a -p as 
follows: 

IMIw-* := sup (g,v), (8) 

v&A 

\\v\\ u s,q=l 

where 1 < q < oo satisfies - + - = 1. We let U~ s,p denote the space U* equipped with 
the U~ s ' p norm. 

Since we can identify elements of U* with elements of U, we can investigate the 
relationship between the U~ 0,p and W 0,p -norms. This will be useful later on in our 
analysis. It turns out that || • \\u-°>p 7^ || • \\u * i n general, but that the following 
equivalence relation holds: 

||w||w-o,p < ||w||wo, P < 2||m|| w -o, p for all u G U. (9) 

To see this, we note that the inequality ||u||w-o,p < ||w||wo, P follows from ([S]) and Holder's 
inequality. To prove the second inequality, we use that fact that, for u EU, 

\\u\\uo, P — sup (u,v) — sup (u,v — v), 

where v — J2f=-N+i v i °f v e ^- 2N ■ Thus, we can estimate 

IM|i/°>p < sup \\v — v\\£9 < 2||li||^-o, P , 

v£R 2N 

\\v\\.q=l 

where we also used the fact that, by Holder's inequality, < \\v\\f? for any v G R 27V . 



Nl»9=l ll"ll»9 =1 



2.5. Projection of non-conservative forces. If we interpret forces as elements of 
U*, then it is natural to consider the following variational formulation of the QCF 
method, 

(^ ci (y) + f,u) = o VueU. (10) 

In other words, (jTOl) requires that J r qc f{y)+f = as a functional in U* . This formulation 
guarantees that the QCF operator has the correct range. 
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To obtain an atom-based description of the equilibrium equations, we explicitly 
compute the representation of Tq^iy) G U* as an element of U (see also (j7j)), that is 
as a zero mean 2iV-periodic vector VuJ r qc f(y), where Vu is defined by 



j=-N+l 

With this notation, the variational equilibrium equations can be understood as pro- 
jected equilibrium equations in atom-based form, 

(V u F q ci(y)) e + fi = 0, -oo < £ < 00. (11) 

The equivalent formulations ffTUl) and ffTTj) define the correct force-based QC method 
for the periodic model problem defined in Section 12.11 

Remark 1. The projection of the QCF equilibrium system is an artifact of the periodic 
boundary conditions. For the displacement boundary conditions that we analyzed 
in [9], or for the mixed boundary conditions that are considered in [6], this projection 
is not necessary. □ 

3. Stability of a uniform deformation 

It is easy to see that, in the absence of external forces, the uniformly deformed lattice 
y = yp is an equilibrium of the atomistic energy as well as the local QC energy, that is 

F a {y F ) = and F c {y F ) = for all F > 0. 

For some values of F, the equilibrium will be stable, by which we mean that the second 
variation 

N 

£"(Vf)[v>,v] = e {0>X + 02VK + ^+i)(^ + ^+i)}> iorueU, 

e=-N+i 

where 

0£:=0"(F) and <^' F := 0"(2F), 

is positive definite, that is, 

£Z(y F )[u,u)>0 WueU\{0}. 

(We note that a second variation, e.g. £"{yp), may be understood either as a bilinear 
form on U or a linear operator from U to U* . It can also be expressed as a Hessian 
matrix with respect to a given basis for the vector space U.) 

In order to avoid having to distinguish several cases, we will assume throughout our 
analysis that F > r*/2, which implies by property (ii) of the interaction potential that 
4>2 F < 0. This assumption holds for most realistic interaction potentials so long as the 
chain is not under extreme compression. 

As above, we can evaluate the second variation of the local QC energy at y = y F , 

N 

£'vi(yF)[u,v}= s Y A Fu' f v'e, 
e=-N+i 
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where A F is the elastic modulus of the continuum model, 

A F :=<p: h {F)=<P" F + A^ F . 

Thus, we say that yp is stable for the local QC approximation if £" cl (y F )[u,u\ > for 
all u G U \ {0}. 

In [8], we have given explicit characterizations for which F the equilibrium yp is 
stable in the atomistic model and in several energy-based QC models. The results for 
the atomistic and the local QC models are summarized in the following proposition. 

Proposition 1 (cf. Prop. 1 and 2 in [8]). Let F > r*/2 then the second variations 
£'l(yF), respectively £ qcl {yp) , are positive definite if and only if 

A F — \ 2 N e 2 (f)2p > 0, respectively if A F > 0, 

where 2 < A a? < n. 

If we denote the critical strains which divide the regions of stability for the atomistic 
and QCL models, respectively, by F* and F*, then a relatively straightforward error 
analysis [8, Sec. 5] shows that F* = F* + 0(e 2 ), that is, the QCL model accurately 
reproduces the onset of a fracture instability. In the following section, we investigate 
whether or not the QCF method has a similar property. 



4. Sharp Stability of the Force-based QC Method 

A trivial consequence of the definition of JF qcf in ([6]) is that y = yp is also a solution 
of the QCF equilibrium equations (fill . 

F^iVF) = for all F > 0. 

(As a matter of fact, this means that the QCF method is consistent; though this is not 
the focus of the present work.) 

To investigate the stability of the QCF method we define the linearized QCF operator 
L qcf ,p:=-F qc{ (y F ):U^ Why 

(L qcitF u,v) := -(^q C f(yF)M,f) for all u, v G U. 

The equilibrium equations for the linearized force-based approximation are then given 
by u G IA satisfying 

(L qc{tF u, v) = (f, v) for all v eU, 

or in functional form 

V u L qci: pu = f. 

We remark that, while L qci)F G L(U, U*), the projected operator VuL qc f t p may be 
interpreted as a map from U to U. 
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4.1. Lack of coercivity. Since the force field J r qc {{y) is non-conservative and the 
linearized QCF operator L qc { p is not the second variation of an energy functional, 
positivity (or coercivity) of L qc f y p may be the incorrect notion of stability for the QCF 
model. Indeed, it turns out that, if N is large, then L qc ^p cannot be positive definite. 



Theorem 2. Let <p" F > and tp'^p ^ 0, then there exist constants C\,C<i which may 
depend on <f>" F and (p'^p, such that, for N sufficiently large and for 2 < K < N/2, 

-dN 1 ' 2 < inf (L qciyF u,u) < -C 2 N 1/2 . 



In [9] , we have shown this result for a Dirichlet boundary value problem. The proof 
carries over from the Dirichlet case almost verbatim and is therefore omitted. As a 
matter of fact, the test function which we explicitly constructed in the proof of Lemma 
4.1 in [9] is already periodic and, after shifting it to have zero mean, can therefore be 
used again to prove Theorem [2J 

Theorem [2] forces us to consider alternative notions of stability. For example, one 
could understand L qc f,F as a linear operator between appropriately chosen discrete 
function spaces, determine for which values of F it is bijective, and estimate the norm 
of its inverse. Physically, this measures the magnitude of the response of the equilib- 
rium configuration to perturbations in external forces, and in Section H~2l we attempt to 
find the largest interval surrounding F — 1 and consider this region to be the approxi- 
mation of the stable region given by operator stability of L qc f F . However, an operator 
can be bijective and its inverse can have bounded norm even when it has negative 
eigenvalues, so for a general equilibrium state such as yp, invertibility alone is not a 
suitable criterion for determining the "physical" stability. To be able to decide whether 
a stable equilibrium of the QCF equations is also stable in a physical sense, we propose 
a notion of dynamical stability in Section 14.41 

4.2. Stability as a linear operator. Since U is a finite-dimensional linear space, the 
choice of topology with which we equip it is unimportant to the question whether L qcijF 
is invertible. However, it has surprising repercussions when we analyze an operator 
norm of the inverse, that is ||£q C f ,p||, in the limit as iV — > oo. 

Our strongest and simplest result is obtained when we view VuL qc { t p as a map from 
U 2 '°° to W - 00 . 

Theorem 3. If \4> F \ — (4 + 2e)|02 F | > 0, then VuL q di,F '-U is bijective and 

1 



2F\ 



Proof. Recalling the definition of L qc f t p, we can rewrite this operator in the form 

^PuL q d,F — <fi F Lx + 4>2f'PuL2, 



10 M. DOBSON, M. LUSKIN, AND C. ORTNER 

where L\ and L 2 are given by 

(L 1 u) g = - e~ 2 (u i+1 - 2u £ + and 
-e~ 2 (u e+2 - 2u e + ui- 2 ), £ ~- 



(L 2 u) 



K,...,K, 



-As (ue+i — 2ug + ue-i), otherwise. 



We note that VuL\ = L\ which is why we have included the projection only in the 
second- neighbor operator. 

The projection of L 2 given by VuL 2 is 



N 



(V u L 2 u) e = (L 2 u) e - - ^ ( L 2 u )j- 

j=-N+l 

We will prove below that 

WPuL 2 \\l(u 2 -°°, u a '°°) < 4 + 2e. 
Assuming that this bound is established, we obtain 

H^w-kqcf.FwIk"" > \<I>f\\\Liu\\£«> - \(j) 2F \\\VuL 2 u\\ ( 
> (|^|-(4 + 2 £ )|^|)|K||, OC) 

which is equivalent to the statement of the theorem. 

To prove (JT21 . we note that, for £ = —K, . . . , K, we have 



(12) 



(L 2 u) e = -K' + i + H + u'U) = -K - K' +1 - 2ul + u'U)- 

Using the first representation of {L 2 u)i above, we immediately see that (for I from the 
continuum region this statement is trivial) 



\{L 2 u)e\ < 4||it"||*» for^=-JV+l, 
From the second representation of (L 2 u)e, we obtain 



N. 



N 



N 



K 



U 



-N+l 



e=-N+i 



e=-K 



U K+1 ' K ' U -K ~ u -K-\i 



and hence, 



\(V u L 2 u)e\ < \(L 2 u) e \ + 



N 



j=-N+l 

< 4||«"||<- + |(K+il + Kl + + Kk\) 

< (4 + 2e)||u"||<«. 

This establishes (fT2l) and thus concludes the proof of the theorem. 



□ 
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Remark 2. With a small modification, Theorem [3] remains true for an arbitrary choice 
of the atomistic region A. The correction 2e then needs to be replaced by n.j£ where 
Hi is the number of interfaces between the atomistic and the continuum region. □ 

Remark 3. Theorem [3] also holds in the case of the artificial Dirichlet boundary 
conditions analyzed in [9]. In that case, the projection Vu is not required and therefore 
the correction 2e does not occur at all. □ 

Theorem [3] is, in many respects, a very satisfactory result. It shows that, except for 
a small error, QCF is stable whenever the atomistic model is. However, the choice of 
function space U 2,0 ° is somewhat unusual, and it is highly unlikely that such a result 
would remain true in higher dimensions, as it requires a regularity that is not normally 
exhibited by linear elliptic systems. 

It is therefore also interesting to analyze the QCF operator as a map from U 1)P to 
U~ l ' p = (U 1,q )*, where 1 < p < oo. However, we saw in [9, Theorem 7.1] for a Dirichlet 
problem that, for 1 < p < oo, the stability of L qc f tF is not uniform in N. The following 
theorem, whose proof is contained in Appendix [Bj establishes the same result for the 
periodic model we consider in the present paper. 

Theorem 4. Suppose that <fi" F > 0, tfi'^p Gl\ {0}, and 1 < p < oo. Then there exists 
a constant C > 0, depending on <j) F and 4>2p, such that, for 2 < K < N — 2, 

H-kqrf.P'lUcw- 1 .*', u 1 ^) > CN 1/p . 



It remains to investigate the case p = oo. The following result is an extension 
of [9, Thm. 5.1] to periodic boundary conditions. Its proof is contained in Appendix IA1 

Theorem 5. If F > r*/2 and <\>" F + 80' 2 V > 0, then 

nr-i II 2 

Theorem [5] establishes operator stability of the L qc f t p operator, uniformly in N, pro- 
vided that 4> F + 802 F > 0. Compared with with Proposition [1] this result predicts a 
significantly smaller stability region than either the atomistic model or the continuum 
model. We employ numerical experiments to see whether the condition <fi" F + 8<p2 F > 
is sharp. 

The norm ^Hl^-i,^ u 1 - 00 ) is difficult to calculate explicitly, so we will estimate 
it in terms of the £°°-operator norm of a related matrix. To that end, we note that, 
according to Lemma [TlJ £ qc f,F can be represented in terms of a conjugate operator, 
-Eqcf.F, by 

(L qc{:F u, v) = (E qc{iF u', v') Vu, v e U. 

An explicit representation of E qc f F is provided in (IT7|) . Formula (|17|) gives an ^ 2Nx2N 
matrix representation for E qcijF such that E qci F e = A F e where e = (1, . . . , 1) T . It thus 
follows that the projected operator VuE qciiF : M. 2N — > M. 2N satisfies 

V u E q A,F ->U, and V u E qc{)F e = 0. 
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Here, and for the remainder of the section, we identify U with the subspace of IR 2Ar of 
zero mean vectors. After these preliminary remarks, we establish the following result. 

Proposition 6. The QCF operator L qc f,F : U — > U* is invertible if and only if 
(VuE^^f + e ® e) G ~$t 2Nx2N ls invertible, and 



§ Halloo < ll-^qcf^lUcW -1 ' 00 , W 1 ' 00 ) — ^ll^l 



where 



T = Vu(V u E qcliF + e ® e)" 1 ^ 



and where \\TWoo denotes the £°° -operator norm ofT. 



Proof. The first statement follows from the discussion above. 

To prove the upper and lower bounds for \\L~^ f pWhiu^^M 1 ' 00 )^ we ^ TS ^ n °t e that, by 
definition of T, it follows that 

TVuE^ F f = P u E qci , F Tf = f for all f eU, 

that is, T = (VuEqc^p)^ 1 on U. In addition, we also have Te = 0. 
Next, we note that 

- — — — = inf sup (L qc{ ^ F v, w) 

\\ L qcf,F\\L(U-^, &><*>) |U/M eW -i „ ™, eM 

IF Wif^- 1 \\w'\\ e i=l 

= inf sup (E qci!F v', w') = — I . 

n„/n e _1 „ ^ w ll^qcf,FlU(w- -°°, 
If llfgo- 1 to' ,i=i 

Since T = ("Pw-Eqci^)" 1 on U, it follows that 

To prove the upper bound, we use fl5]) to estimate 

M lk(W- '°°, M°>°°) = SU P "TTTTl — SUp j-r— n < 2 SUp = 2||J ||oo. 

... II >- II . .... 2JV 



To prove the lower bound, we first note that TVu = T. We will also use the fact that 
H^w/ll^f < 2||/||^«> for all / G R 2N . Employing also (jHD again, we can deduce that 

|m , \\TV u f\\u°.°° . IIT/I^go 
||r|| z(w -o,oc, w o lOC) = sup — — > sup = sup ° 

/gr 2JV Ikw/IIw- - 00 /gm 2JV Ikw/IIw ' 00 /gr 2JV Ilnw/lki° 
\\Tf\Ur 1 IgTfe i im 



jm 2N Z WJ \\t? z /gm 2JV 

The penultimate equality holds because Vuf = implies T f = 0. □ 

In Proposition [6] we have reduced the estimation of the operator norm of L~ { F to 
the computation of the £°°-operator norm (which is simply the largest row sum) of a 
matrix T G R 2Nx2N , which is explicitly available (note that V u = I-|e®e G R 2iVx2Ar ). 
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Stability Bound for L qc f p 

10 3 



10 2 

J 

10 1 



0.2 0.4 0.6 0.8 1 

1 + 40' 2 ' f /0f 

Figure 1. Computation of HTH^ where T = Pw(Pt/£qcf,F + e<S>e)~ 1 'P w , 
which gives lower and upper bounds for \\L~^ pW^u- 1 ^. (cf- Propo- 
sition [6]). The graphs indicate that L qc ^p is stable as an operator from 
U 1 ' 00 to W^ 1 ' 00 , uniformly in N, for all macroscopic strains F up to the 
critical strain for QCL and QNL. 

In Figure [TJ we plot the norm of T as a function of Ap/(j) F = 1 + A^p/^p. We 
clearly observe that L qc f,F is in fact stable for all macroscopic gradients F for which 
Ap > 0, that is, the bound required in Theorem[5]is not sharp. Moreover, the numerical 
experiments shown in Figure [TJ support the following conjecture. 

Conjecture 7. If (fi F + A^p > 0, then 

where r\ does not depend on N or K, but r/^l + 4^f j ^00 as 1 + 4^f — > 0. 

In fact, the numerical experiments suggest that ||£q C f f ||l(w-i,°°, m 1 * 00 ) grows faster 
than ,„}.,,, , which would imply that an estimate such as the one in Theorem but 
with the constant 8 replaced by 4 would be false. 

4.3. The quasi-nonlocal coupling method. In preparation for the following sec- 
tion, where we introduce another notion of stability for the QCF method, we review 
a popular energy-based coupling method. In the next section, we will make numerical 
comparisons between this method and the QCF method. 

The quasi-nonlocal quasicontinuum approximation (QNL) [26] was derived as a mod- 
ification of the energy-based QC approximation [20] in order to correct the incon- 
sistency at the atomistic-to-continuum interface [5,24]. In the case of next-nearest 
neighbour pair interaction, the QNL method can be formulated as follows. Nearest 
neighbor interaction terms are left unchanged. A next-nearest neighbor interaction 
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term <f>(e 1 (ye+i — Ve-i)) is left unchanged if atom £ belongs to the atomistic region, 
but is replaced by a Cauchy-Born approximation 

</>{e- 1 (y l+1 - y<_i)) « \ [<j>{2y' t ) + 0(2^ +1 )], if I e C. 

This process yields the QNL energy functional 

N 

We remark that the QNL method is consistent for our next-nearest neighbour pair 
interaction model, and in particular, y F is an equilibrium of the QNL energy functional 
in the absence of external forces. Moreover, in [8] we have established the following 
sharp stability result for the QNL method, which shows that the QNL method is 
predictive up to the limit load for fracture. 

Proposition 8 (Proposition 3 in [8]). Suppose that F > r*/2 and that K < N — 1, 
then £qni(zAF) is positive definite in U if and only if A F > 0. 

4.4. Dynamical Stability. We have pointed out in Section H~T1 that operator stability 
for L qcf >F cannot guarantee that the equilibrium y F is a stable equilibrium of the 
atomistic model (e.g., a local minimum). To obtain at least a theoretical methodology 
to determine stability of yp from the QCF operator alone, we propose the notion of 
dynamical stability. The dynamical system 

u(t) + V u L qci , F u(t) = 0, 
u(0) = m , «'(0) = 0, 

has a unique solution u G C°°([0, +oo); IX). We call this dynamical system stable if 
there exists a constant C, independent of N, such that 

IKt)||*§ < C\\u \\p e Vt > 0, Vm e U. (13) 

This condition can be best understood in terms of the spectrum of P«£ qc f,F- In numer- 
ical experiments, which are shown in Table [H we have made the surprising observation 
that VuLqct,F and S'^iyp) appear to have the same spectrum. This has led us to make 
the following conjecture. 

Conjecture 9. For all N > 4, 1 < K < N, and F > 0, the operator VuL qc f, F is 
diagonalizable and its spectrum is identical to the spectrum of E'LJ^yp). 

Since S^ nl (y F ) is positive if and only if A F > (cf. Proposition [8]), the validity of the 
conjecture would imply that L qcijF has positive real eigenvalues if and only if A F > 0. 

To see how this observation implies dynamical stability ([TBI , let V denote the 
matrix whose columns are the eigenvectors for VuL qc f tF . Then V has full rank and 
V^VuLqcf^V is a diagonal matrix with the eigenvalues of L qc f tF on its diagonal. If 
we define z{t) = V~ l u{t), then 

z(t) + V- x V u L qci ,pVz(i) = 0, 
z(0) = V^uo, z'(0) = 0. 



SHARP STABILITY ESTIMATES FOR QCF METHODS 



15 



N 







-0.1 



-0.15 



-0.2 



-0.25 



2F — 



50 
100 
150 
200 
250 
300 












1.19e-010 
6.97e-010 
2.05e-009 
4.44e-009 
8.25e-009 
1.62e-008 



9.93e-011 
6.19e-010 
1.83e-009 
3.12e-009 
6.38e-009 
1.15e-008 



7.31e-011 
4.71e-010 
1.31e-009 
2.90e-009 
6.38e-009 
9.98e-009 



6.64e-011 
3.16e-010 
1.23e-009 
2.10e-009 
3.96e-009 
8.86e-009 



Table 1. The spectra of Pw£ qc f,F and S'^ nl (yp) are computed for in- 
creasing N, for K = N/2, for <p F = 1, and for different values of (p'^p. 
The table displays the £ 2 norm (not scaled by e) of the ordered vectors 
of eigenvalues. The column for = is identically zero since, in this 
case, the two operators coincide. All other entries are zero to numerical 
precision of the eigenvalue solver. 



N 


<&F = 


-0.1 


-0.15 


-0.2 


-0.24 


10 


1.00 


1.4563 


1.7607 


2.3770 


5.4398 


30 


1.00 


1.5242 


1.9441 


2.8049 


6.2876 


90 


1.00 


1.5794 


2.0537 


3.0726 


7.4324 


270 


1.00 


1.6049 


2.1095 


3.1510 


8.2878 


810 


1.00 


1.6136 


2.1191 


3.2021 


8.3863 


2430 


1.00 


1.6139 


3.7477 


3.2075 


8.5968 



Table 2. The table shows the condition number cond(V) for increasing 
values of N, for K = N/2, for (f> F = 1, and for different values of 2 ' F . In 
each column, the computed condition numbers appear to approach an 
upper bound. 



where the condition number of V is defined as usual by cond(V) = \\V\\ || V _1 ||. Hence, 
we see that, subject to the validity of Conjecture [HI the L qc{ F operator satisfies (TT31 
with constant C = cond(V). To make this stability independent of N, we require that 
cond(K) is bounded as iV — > oo. This is the subject of further numerical experiments 
displayed in Table [2j They suggest that this is indeed true if and only if Ap > 0. 

Conjecture 10. Let V denote the matrix of eigenvectors for the force-based QC 
operator VuLqcf^F- If Ap > 0, then cond(V) is uniformly bounded in N. 

Conjectures H2 and [TDl supported by the results of the numerical experiments that 
we have presented in Tables [TJ and [2], imply that VuL qc f^F is indeed dynamically stable 
for Ap > 0, with a stability constant that is uniform in N. 
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Conclusion 

We propose that a sharp stability analysis of atomistic-to-continuum coupling meth- 
ods is an essential ingredient for the evaluation of their predictive capability, as impor- 
tant as a sharp consistency analysis. In the present paper, we have established such a 
sharp stability analysis for the force-based QC method, for a simple one-dimensional 
model problem. We have analyzed three notions of stability: 

(i) Positivity (coercivity) is generically not satisfied. 

(ii) Operator stability, uniformly in the size of the atomistic system, holds only with 
an appropriate choice of function spaces. It does not hold for several natural 
choices. 

(iii) Dynamical stability is satisfied up to the critical load. This result is based on 
the numerical observation that the spectra of the QCF and QNL operators 
coincide. 

Positivity and dynamical stability are equivalent for energy-based methods, and 
under suitable conditions and choices of function spaces they imply operator stability. 
However, the fact that the QCF method is non-conservative and gives rise to non- 
normal operators, leads to a much richer mathematical structure. 

Finally, we stress once again that, while the QCF method is possibly the simplest 
force-based multi-physics coupling scheme, we believe that similar observations can be 
made for other force-based hybrid methods, such as FeAt [15], CADD [25] or brutal 
force mixing [2]. 

Appendix A. Proof of Theorem O Stability of L qcf F 

Theorem [5] states that, if 4>" F + 80 2 ' F > 0' then L qci F is stable as an operator from 
U 1 ' 00 to W -1 ' 00 , uniformly in N. 

The proof of this statement uses a variational representation for the QCF operator, 
which we derived in [9], and which is also valid for periodic boundary conditions: 

where the three operators L 1; L™ s , : U —>■ U* are given by 
(Ltu,v) = (u',v'), 

-K K N 

(L T 2 eg u, v) = eY^ ^ u We + £ Yl ( u <-i + 2u 't + u e+^ v t + £ Y1 Au '^ 
e=-N+x e=-K+i e=K+i 

(L S 2 n& U, V) = (u'_ K +l ~ 2u '-K + u '-K~l) v ~K - iu' K+2 - 2U K+1 + U' K )V K . 

We omit the proof of this representation which is a straightforward summation by parts 
argument and carries over verbatim from [9]. Upon defining 

cj)" F u' t + 02F«-i + 2u' e + u' i+1 ), £ = —K + 1,...,K, 



a e (u') 



b'p + 402^)^) otherwise, 



as well as 

a K {u ) = (p 2F { u K+2 ~ ' Zu k+i + u k), 



- 2u'tr, n + u'tA, and 



a- K (u') = (p2 F (u'_ K+1 - 2u'_ K + u'_ K _ l ), ^ ^ 
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we can rewrite this representation as 

(L qcf!F u,v) = (a(u'),v') + ee- K (u')v- K - a K (u')v K - 
Using the periodic heaviside function h EU given by 

ki -\ -1(1 + ^) - |, l<0, W 

and setting he = he-i, the point evaluation functional «ho 0) «gW, can be represented 
by 

v o = (h',v) = — (h, v') for all v E U. 
Combining these observations, we obtain the following result. 

Lemma 11. L qc f^ F can be written as 

(L qcl:F u, v) = (E qci:F u', v') for all u,v EU, (16) 

where 

Eq C f,FU e = &e(u') - a_ K {u')hi +K _i + ct K {v!)h,>_ K _ x , (17) 
for a, h, and a±K as defined above. 

Even though the variational representations of the Dirichlet case and the periodic 
case are the same, we cannot translate the proof for inf-sup stability that we used 
in [9], as it required a matrix representation that is unavailable for periodic boundary 
conditions. Instead, we will compute a fairly explicit characterization of L~^ f F to 
estimate its norm directly. It is most convenient to do so if we define an equivalent 
norm on U~ l,oc . Note that L\ : U — > U* is bijective, and hence we can define 

Ibllw-i.oc = ll-^r^llw 1 ' 00 ioigEU*. 

Lemma 12. For all g EU* , it holds that 

1 

2 IMIw- 1 . 00 — IMIw- 1 ' 00 < IMIw- 1 ' 00 - 
Proof. Let z = L± x g, that is, 

(z',v') = (g,v) VvEU. 
Taking the supremum over v with \\v'\\ei = 1 we obtain the second inequality 

llgllw-!. 00 < Ik'lkf — Ibllw- 1 - 00 

by Holder's inequality. 

The first inequality follows from the fact, which is proved below, that 



i ii „/ 



z'Weoo < sup (z',v') \/z E U. (18) 
IKIU=i 



Namely, this implies that 



iIMIw-l- = \\\ z 'hr ^ SU P ( z ', v ') = SU P (g,v) = \\g\\u- 1 - 
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To prove (fT8lh we fix z G U and let £1,^2 be such that z' £i = \\z'\\£&> and z[ 2 < 0. 
(A similar argument can be used if z' £i = — \\z'\\£oa.) We obtain ( TTHl) from the fact that 



7T \ \ Z foe < (z', v') where v G IA is defined by 



v' 



i if£ = 4, 

if* = 4, □ 
otherwise. 



Corollary 13. Suppose that F is such that L qc f,F :U — > W* invertible, then 

\\(L 1 1 L qc f )F ) 1 \\l(U 1 >°°, U 1 ' 00 ) < H-^qcf.FlUcW -1 ' 00 , 

— 2|| (L 1 1 L qc f,F) 

Proof. Using Lemma [T2] twice, we can prove the following bound, 

M/ r -l r TZTn = in £ 11-^1 1 ^qcf,F«||w 1 .°° = in 7 f ||-^qcf,F^||^-i, 

II (A ^qcf,Fj |U(Wl.°°, „ n" 6 " 1 1 |, w6M , 

11^11^1,00 — J- 11^11^1,00 — -L 

> inf ||L^ 1 Lq C f j Fw||wi,oo 



" , . II 1 ^qcr.j- iir-1 
i&A \\Lj c 

iiwi 

which gives the first stated inequality. The second inequality follows from a similar 
argument. □ 

Corollary [131 shows that we can bound the operator norm \\L~ t f ||l(w-i,°°, u 1 ' 00 ) m 
terms of ||(Z/J" ^qcf,^) - The latter operator norm can be computed using 
the formula 

\\(L 1 1 L qc{:F )~ 1 \\uu 1 ' o °, u 1 ^) = \ inf -^qcf,Fw)'||i2° \ ■ (19) 

lu'lUoo =1 

II lltg 

In the next lemma, we establish an explicit representation of L^ 1 L qc ip which will 
subsequently allow us to construct upper and lower bounds for (fl9|) . 

Lemma 14. Let z = L^ l L qc {^pu, then 

H = W ) - -j<?2F\ u -K ~ U -K+1 -U K + U K+1 \ 

- a- K {u')hz +K _i + a K {u')h f ^ K ^i, 
where a, h, and a±K are defined above. 

Remark 4- We note that the term \e{u'_ K — u'_ K+l — u' K + u' K+1 } is the average of 
a, and the function h is a periodic heaviside function defined in (1151) . □ 



Proof of Lemma [X|j The function z is the solution of the variational principle 

(z',v') = (L qcftF u,v) = (E qciiF u',v') 
where E qci F is defined in ( TTT1) . and is given by 

E qc{iF u' e = o- £ (u') - a- K (u')h e+K _i + a K {u')ht- K -x- 
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We note that a function w G R is a gradient, that is, w = v' for some v G U, 
if and only if Yli=-N+i w t = ®- Hence, we obtain z' = E qc { tF u' — E qc f tF u' where 



Eqci,Fu' := ^J2e=-N+i Eqci,Fu' t Since h has zero mean, we only need to compute a, 
1 N A F N (f)" K 



a ' = 2iV ^ "* 2iV ^— ' ~* ' 2iV 

£=-JV+l £=-AT+i £ = _jf + i 



Since u is periodic, u' has zero mean, and hence the first sum on the right-hand side 
vanishes. The second sum has telescope structure, and we obtain 



e 



E qc{)F u' = a = -4>2 F {u'_ K - u'_ K+1 - u' K + u' K+l ). 
This concludes the proof of the lemma. □ 

We are now ready to conclude the proof of Theorem [5j 
Proof of Theorem^ We set z = L qc { tF u and use Lemma fl4l to deduce the bound 
\W\\t? > \W(u')\U°° - 2£|(/)2 F |||« / ||£ ? c 

- max(\a- K (u')\, \a K (u')\)max(\he +K -i\ + \he- K -i\)- ^ 
To bound the first term on the right-hand side, we note that 

We(u')\ > 4> F \u' e \ +4(f)2 F \\u'\\e r , 

which immediately implies 

\\(j(u')\\ er >A F \\u'\\ er . (21) 
To bound the third term on the right-hand side of (TSUI) , we crudely estimate 

max + \h t -K-i\) < 1 - h, 

e=~N+i,...,N 2 

which is true whenever K > 1, and deduce from (|14p that 

max(\a- K (u')\, \a K (u')\) < 4|0' 2 ' F |. 

The additional term — |e cancels with the second term on the right-hand side of (120]) . 
so that we obtain 

\\z'\\e r > (4>' F + 8(f)2 F )\\u'\\ i o . 
Employing Corollary [13] and Formula (119j) . we obtain Theorem [51 □ 

Appendix B. Proof of Theorem m Instability of L qcitF 

We now prove Theorem H] on the instability of L qc £ tF as an operator acting between 
U x ' p and U~ x *, 1 < p < oo. The bound \\L^ F \\ L{ u-l P)A x,v) > CN l / p follows from the 
following lemma. 

Lemma 15. Suppose that <fi F > 0, tfi'^p G R \ {0}, and p, q G R satisfy 1 < p < oo, 
1 < o < oo, and - + - — 1. TTien t/iere exists a constant C > stzc/i i/iai 

inf sup (L qcf)F t», tu) < CiV^ 1/p . 
ll u lu? =1 IKIU=i 
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Proof. We recall from Lemma [H] that we can represent L^pv in the form 

(L qcf)F v , w) = (E qc { tF v', w') Vu> e U, 



(22) 



where 



and where 



E qc {,Fv'i = ae(v') - a_ x (i/)/i £+x _i + a K (v ')h t -K-i, 



aAv 



a K {v') = 4>2 F (v' K+2 - 2v' K+1 + v' K ), 

a- K (v') = 4>2 F ( V '-K+1 - 2v '-K + v '-K-l), 

>o, 

< 0. 



£=-K + l, 
otherwise, 



,K, 



he 



\{l-e£)-\ 
-Ul + et) - 6 



2\ 1 4' " 

We choose v G U with derivative given by 

r o, 



6< F ' 
A F 

A F 



hp 



-K-li 



e = K-i, 

£ = K, 

e = K + i, 

£ = K + 2, 
otherwise. 



Such a representation is possible if and only if the vector (v' e )^L_ N+1 defined above has 
zero mean. To see that this holds, we use the symmetry of he to calculate 



N 

E «; 

e=-N+i 



he-i<-i 

l^K-l,K,K+l,K+2 



0. 



If we insert v into the equations above, we find that 

a- K (v') = 0, a K (v') = —Ap, 



and 



aAv 



{(j) F + 2$ F )/l_ 3 + $A-4, 



b' F A F 



4 2 



4 2 



Aphp 



-K-1-, 



£ = K-2, 
£ = K-1, 

£ = K, 

£ = K+1, 

£ = K + 2, 
otherwise, 
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which implies that 



E qc i,FVe 



-20^-3 + <% F h-A, 

% F h- 3 - \A F - A F h_ 2 , 

A 2 



A F h , 



0. 



e = K -2, 

£ = K - 1, 

i = K, 

£ = K + 1, 

£ = K + 2, 
otherwise. 



Note that all the terms above are bounded in absolute value, independently of N and 
K. 

Inserting these formulas into (1221) . applying Holder's inequality, and using the fact 
that E qc f^pv[ is nonzero for only five indices, we obtain 

(L qc{>F v, w) = (E qc{>F v', w') 

< \\Eqci,FV \\fP 

jl/p 

' J qci,l """" 



< e 



i/p 



5 H-E'qcf.F'w'IlL 



\w u 



< Ce 1/p 



w 



It remains to show that Hf'll^ is bounded below as iV — >• oo. As a matter of fact, it 
can be seen from the definition of v' s that 



which gives 



\v,\ > 



> 



for j = K + 1- N/2,...,K -1, 



=K+l-N/2 

Thus, replacing v by vj \v'\g> gives the desired result. 



K-l 



i Vp 



\ [(N/2 - 2)e] 



i/p 



□ 
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